library(SDMTools)
wd='/home/jc148322/rob/';setwd(wd)

#read in the asciis
bdall=read.asc('Bdall.asc')
bdwet=read.asc('Bdwet.asc')
bc12=read.asc('bc12.asc')
bc01=read.asc('bc01.asc')
veg=read.asc('veg.asc')


######################################
#create the data
pos = as.data.frame(which(is.finite(bdall), arr.ind=TRUE))
pos$lat = getXYcoords(bdall)$y[pos$col]; pos$lon =
getXYcoords(bdall)$x[pos$row] #append the lat lon
pos$bdall= extract.data(cbind(pos$lon,pos$lat),bdall)   #extract the data
pos$bdwet= extract.data(cbind(pos$lon,pos$lat),bdwet)
pos$bc01= extract.data(cbind(pos$lon,pos$lat),bc01)
pos$bc12 = extract.data(cbind(pos$lon,pos$lat),bc12)
pos$veg = extract.data(cbind(pos$lon,pos$lat),veg)


#read in the threshold to create sigdiff
threshold=read.csv('maxentResults.csv')
threshold=threshold$Minimum.training.presence.logistic.threshold
bdall[which(bdall<threshold[1])]=0; #set values below threshold to 0
bdwet[which(bdwet<threshold[2])]=0; #set values below threshold to 0
sig.diff=SigDiff(bdall,bdwet)
pos$sig.diff=extract.data(cbind(pos$lon,pos$lat),sig.diff)
write.csv(pos,'coord.csv')

#read in the frog data and subset it to columns of interest
pos=read.csv('coord.csv',as.is=TRUE)
frogdata=read.csv('torrent.csv', as.is=TRUE)
frogdata=frogdata[,c('Species','bc01','bc12')]
lorica=frogdata[which(frogdata[,1]=='LITLORI'),]; lorica[,2]=lorica[,2]/10
nannotis=frogdata[which(frogdata[,1]=='LITNANN'),]; nannotis[,2]=nannotis[,2]/10



######################################
#Create the plot
#edit and simplify the data
tdata = cbind(pos$bc01,pos$bc12); #
tdata[,1] = round(tdata[,1],1); tdata[,2] = round(tdata[,2]/25)*25
#round bc01 and bc12
tdata = unique(tdata) #find only unique combinations of rounded bc01 and bc12

pos$bdall[which(pos$bdall<threshold[1])]=0; #set values below threshold to 0
pos$bdwet[which(pos$bdwet<threshold[2])]=0;
# pos$bdall[which(pos$bdall>0)]=1 #set values above threshold to 1
# pos$bdwet[which(pos$bdwet>0)]=1
# tdata.all = unique(cbind(pos$bc01[which(pos$bdall==1)],pos$bc12[which(pos$bdall==1)]))#find only unique combinations
# tdata.wet = unique(cbind(pos$bc01[which(pos$bdwet==1)],pos$bc12[which(pos$bdwet==1)]))#find only unique combinations

tmp=pos[,c('bc01','bc12','bdwet')]
tmp=tmp[which(tmp$bdwet>0),]
tmp$bdwet=round(tmp$bdwet*100)


tmp2=pos[,c('bc01','bc12','bdall')]
tmp2=tmp2[which(tmp2$bdall>0),]
tmp2$bdall=round(tmp2$bdall*100)


#make the plot
# cols=colorRampPalette(c("tan","forestgreen"))(5)
palette(colorRampPalette(c("tan","forestgreen"))(100))
png('compare_frogs_test.png', width=16,height=8,units='cm', res=300,pointsize=5)
	par(mfrow=c(1,2), mar=c(5,5,1,1))
        plot(tdata[,1],tdata[,2], xlab='Annual mean temperature',ylab='Annual rainfall',xlim=range(tdata[,1],na.rm=T),ylim=range(tdata[,2],na.rm=T), type='n',cex.lab=1, cex.axis=1, font.lab='2')
        #points(tdata[,1],tdata[,2], col='grey', pch=15)

        #points(tdata.all,col='tan', pch=15)

        points(pos$bc01[which(pos$bdwet>0)],pos$bc12[which(pos$bdwet>0)],col=tmp$bdwet,pch=15)

        # points(nannotis[,2],nannotis[,3],col='red',pch=15)
        # points(lorica[,2],lorica[,3],col='blue',pch=15)
		text(22,6000, 'bdwet')
		
		plot(tdata[,1],tdata[,2], xlab='Annual mean temperature',ylab='Annual rainfall',xlim=range(tdata[,1],na.rm=T),ylim=range(tdata[,2],na.rm=T), type='n',cex.lab=1, cex.axis=1, font.lab='2')
        #points(tdata[,1],tdata[,2], col='grey', pch=15)

        #points(tdata.all,col='tan', pch=15)

        points(pos$bc01[which(pos$bdall>0)],pos$bc12[which(pos$bdall>0)],col=tmp2$bdall,pch=15)

        # points(nannotis[,2],nannotis[,3],col='red',pch=15)
        # points(lorica[,2],lorica[,3],col='blue',pch=15)
		text(22,6000,'bdall')
dev.off()

######################################
#make the image

out.range=range(sig.diff, na.rm=T)


Colormap=c("grey",colorRampPalette(c("tan","forestgreen"))(100))
#Colormap2=c(colorRampPalette(c("tan",'burlywood'))(length(unique(ctan))),colorRampPalette(c("olivedrab2","olivedrab"))(length(unique(cgrey))),colorRampPalette(c("forestgreen",'darkgreen'))(length(unique(cgreen))))
cols=c('forestgreen','grey','tan')
pnts=cbind(x=c(429509,446345,446345,429509),y=c(8205059,8205059,8128594,8128594))

png ('sigdif_bd_testing.png',width=3*dim(bdall)[1],
height=1*dim(bdall)[2], bg= "white", pointsize=100)
    par(mfrow=c(1,3))

        # image(bdall, col=Colormap, ,xlab="", ylab="")
        # legend.gradient(pnts,col=Colormap, c(0,1), title='Original')
        # image(bdwet, col=Colormap, ,xlab="", ylab="")
        # legend.gradient(pnts,col=Colormap, c(0,1), title='New')

        image(sig.diff,col=cols, zlim=range(sig.diff,na.rm=T),xlab="", ylab="")
        legend.gradient(pnts,col=cols,round(range(sig.diff,na.rm=T),digits=1), title='SIG.DIFF')
		
		ImageDiff(sig.diff,sig.levels=c(0.025,0.975),tcol=cols)
		
dev.off()







